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^ Abstract 

\o 

O We present a detailed Monte Carlo study coupled to a likelihood analysis of the po- 
^"^ tential of next generation ultra-high energy cosmic ray experiments to reconstruct 
properties of the sources and the extra-galactic magnetic field. Such characteristics 
0"\ are encoded in the distributions of arrival time, direction, and energy of clusters 
"f-j of charged cosmic rays above a few 10 19 eV. The parameters we consider for recon- 
Q-6truction are the emission timescale, total fluence (or power), injection spectrum, 
O and distance of the source, as well as the r.m.s. field strength, power spectrum, and 
£j coherence length of the magnetic field. We discuss five generic situations which can 
^ be identified relatively easily and allow a reasonable reconstruction of at least part 
^ of these parameters. Our numerical code is set up such that it can easily be applied 
• ■'-J to the data from future experiments. 

^5 PACS numbers: 98.70.Sa, 98.62.En 
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1 Introduction 



The origin and the nature of ultra-high energy cosmic rays (UHE CRs), with energies E > 10 EeV (1 EeV = 
10 18 eV), are still unknown despite several generations of experiments, most notably the Haverah Park [1], 
the Akeno Giant Air Shower Array (AGASA) [2,3], and the Fly's Eye [4] experiments. Data from this 
latter seem to indicate that the UHE CR component is mainly composed of protons [4]. At these energies, 
protons cannot be confined within the Galactic magnetic field. Thus, the isotropy of the arrival directions 
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of most of the observed UHE CRs [3] , or at least the absence of a significant correlation with the Galactic 
plane [5], suggests that UHE CRs are extra-galactic in origin. 

However, such protons would leave a distinct signature in the energy spectrum, the so-called Greisen- 
Zatsepin-Kuzmin high energy cut-off [6] (hereafter GZK cut-off) around E ~ 70EeV, due to pion produc- 
tion on the cosmic microwave background by nucleons with E > 70 EeV. There is no strong experimental 
evidence for this cut-off and the detection of particles with energies as high as E ~ 300 EeV cannot be 
easily explained in this frame. No compeling astrophysical candidate for the source of the highest energy 
events could be found within ~ 100 Mpc [7,8], although photopion production limits the range of nucleons 
with E ~ 100 EeV to about 30 Mpc. Heavy nuclei would be disintegrated over similar or slightly larger 
distances [9], and similar problems arise for the less likely option of a 7-ray [8], for which the effective 
attenuation length in electromagnetic cascades lies between 1 and 20 Mpc, depending on the poorly known 
strength of the universal radio background. Finally, neutrino primaries in general imply too large a flux 
because of their small interaction probability in the atmosphere [10]. 

As to the theoretical models of the origin of UHE CRs, the most conventional scenario involves first-order 
Fermi acceleration of protons in powerful astrophysical shocks, for instance in the hot spots of radio- 
galaxies [11]. More recently, it was suggested that protons could be accelerated up to E ~ 10 21 eV in 
fireball models of cosmological 7— ray bursts [12-15]. In order to reconcile the observed rates of UHE CRs 
and cosmological 7— ray bursts within D ~ 30 Mpc, however, the arrival time of UHE CRs would have 
to be spread over At > 50yrs, for instance through deflection in large-scale magnetic fields [12,13,16]. As 
another class of models, topological defects, possible relics of early Universe phase transitions, could release 
supermassive "X" particles with mass around the Grand Unification Scale, through physical processes such 
as collapse or annihilation [17]. These X particles would subsequently decay to jets of UHE CRs, with a 
likely dominance of 7— rays above ~ 50 EeV, and an energy spectrum significantly harder than in the case 
of shock acceleration [18,19]. 

Charged UHE CRs, such as protons or electromagnetic cascades initiated by a 7— ray primary, are subject 
to energy-dependent deflection, and hence energy-dependent time delay, in large-scale magnetic fields. The 
r.m.s. strength B rms and the coherence length l c of extra-galactic magnetic fields are thoroughly unknown, 
although they are bound by Faraday rotation data to B rms l^ 2 < 10 _9 GMpc 1 ^ 2 [20]. Different authors 
have proposed to use UHE CRs to probe extra-galactic magnetic fields; some rely on the magnitude of the 
time delay and the deflection [21,22], or on some features of the angle-time-energy images of UHE CRs 
[23,24], or even on synchrotron loss signatures in the energy spectrum of electromagnetic cascades [25]. In 
a previous paper, we discussed how information on both the extra-galactic magnetic field and the origin of 
UHE CRs could be left in angle-time-energy images of clusters of proton UHE CRs [24]. We applied this 
study to a maximum likelihood analysis of the three pairs of UHE CRs [26], that were reported by the 
AG AS A experiment [3] . In order to do so, we devised a Monte-Carlo code that follows the propagation of 
UHE protons and calculates a likelihood as a function of the parameters characterizing the origin of these 
UHE CRs and the intervening magnetic fields. 

Future large scale experiments [27], such as the High Resolution Fly's Eye [28], the Telescope Array [29], 
and most notably the Pierre Auger project [30], should allow to detect clusters of > 20, and possibly 
more, UHE CRs per source, if the clustering suggested by the AGASA results [3] is real. The recently 
proposed satellite observatory concept for an Orbital Wide-angle Collector (OWL) [31] might even allow 
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to detect clusters of hundreds of events by watching the Earth's atmosphere from space. In the present 
paper, we wish to examine how magnetic fields could affect the observations of clusters of UHE CRs 
by future large-scale experiments. We thus assume that UHE CRs are indeed dominantly protons, and 
that the magnetic fields are strong enough to influence their propagation (see below). We then simulate 
the injection, the propagation, and the detection of UHE CRs originating from a given source. Finally, 
we perform a maximum likelihood analysis on these clusters of typically 20 — 50 particles, and attempt 
to reconstruct the physical parameters describing the source and the magnetic fields. We describe the 
simulations in Section 2, and discuss the reconstruction of the different parameters in Section 3; we briefly 
summarize our results in Section 4. We use natural units, h = c = 1, throughout the paper. 



2 Simulation of clusters of Ultra-High Energy Cosmic Rays 

Protons of ultra-high energy are subject to the following physical processes: energy loss through pair 
production and photo-pion production (the latter for E > 70 EeV) on the cosmic microwave background, 
and deflection in the extra-galactic magnetic field. In photo-pion production, a proton may be converted 
to a neutron, that either turns back into proton through photo-pion production, or decays to a proton on a 
distance ~ 1(E/10 20 eV) Mpc for a neutron energy E. Pair production is treated as a continuous energy loss 
[32]. We treat photo-pion production as a stochastic energy loss; it is important to do so, as the stochastic 
nature of this process imprints significant scatter in arrival time and energy for UHE CRs above the GZK 
cut-off, as discussed in Ref . [26] . We model the extra-galactic magnetic field as a gaussian random field, with 
zero mean, and a power spectrum given by (B 2 (k)) oc k nB for k < 2ir/l c , and (B 2 (k)) = otherwise. The 
cut-off, l c , characterizes the coherence length of the field. The field is actually calculated on a grid of inter- 
cell separation a and is tri-linearly interpolated between the lattice points such that l c ~ a/n effectively. 
The amplitude of the field is normalized to the r.m.s. strength B 2 ms = V/(2n) 3 J d 3 k£? 2 (k), and our model 
for the extra-galactic magnetic field is thus described by the three parameters ris, l c , and B TUIS . Fiducial 
values for these parameters are % ~ 0, l c > 100 kpc, and -B rms ^ 10~ 9 G. This statistical description of 
the field allows to treat deflection of UHE CRs in the most general case, as discussed in Ref. [26]. We also 
note that any relative motion between observer and source with relative velocity v would introduce effects 
only on timescales larger than l c j v which is much larger than delay times and experimental lifetimes. It is, 
therefore, justified to assume a stationary situation. 

The numerical code that we use to follow the propagation of UHE protons in an extra-galactic magnetic 
field is described in detail in Ref. [26]; here, we summarize its main features. Protons are injected with a flat 
energy spectrum, and propagated in a given direction in the extra-galactic magnetic field over a distance 
D, from the source to the detector. Due to the stochastic deflection, care has to be taken in how one states 
whether different UHE CRs, that have followed different paths, actually reached the same detector, or not 
[26]. During their propagation, UHE CRs acquire energy-dependent deflection 6 E and time delay te- With 
a given sample of nucleons, one can construct different histograms, in time, angle, and energy, for different 
values of the differential injection index 7, and of the fluence N . Histograms are then smeared out in energy 
with AE/E = 0.14, i.e. a high resolution typical of future large-scale UHE CR experiments; histograms 
are also convolved in time with a top-hat of width T$, in order to simulate emission of particles at the 
source over a timescale T$. Once the histogram is obtained for different values of the above parameters, 
clusters of UHE CRs can be obtained by picking at random a time window of length T obs ~ 5 yr, which 
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corresponds to the lifetime of the experiment, and dialing Poisson statistics over the histogram. We do so 
in order to simulate UHE CR clusters of events. 

Conversely, one can use the above code to perform Monte-Carlo simulations of UHE CR injection, prop- 
agation, and detection, and calculating a likelihood of a given histogram for a given cluster of events, 
where the histogram, hence the likelihood, is a function of the physical parameters described above. The 
likelihood is calculated in the standard way for each observed event cluster, using Poisson statistics, 

C(r W o,T s ,D n ,N ,l c ,n B ) = ( 1] e " Pj ^hr) > (!) 

\j=l,jV n U)-/ 



where pj is the predicted number of events in cell j, and n(j) is the number of observed events in cell j for 
the cluster under consideration. Each cell is defined by a time coordinate and an energy. The time-energy 
histogram is binned to logarithmic energy bins of size 0.05 in the logarithm to base 10 (as opposed to 0.1 
in Ref. [26] to account for improved energy resolution of future experiments), and to 0.1 yr in linear time 
bins. The product in Eq. (1) extends over all energy bins (from 10 1 ' 5 EeV to 10 4 EeV) and over all time 
bins within an observational time window of length T obs ; we took T obs ~ 5yr as a projected lifetime of a 
next generation experiment such as the Pierre Auger Project. The brackets in Eq. (1) indicate that the 
likelihood has already been averaged with equal weights over the position of the observational time window 
on the time delay histogram of the UHE CRs, as well as over different realizations of the extra-galactic 
magnetic field between the source and the observer. 

The next step is to attempt to reconstruct the parameters in Eq. (1) from the maximum of the likelihood. 
Future experiments are expected to produce as many as > 100 particles with E > 50 EeV, if the AGASA 
pairs are real. In the present work, we prefer to remain conservative, and we simulate clusters of 20 — 50 
particles with E > 30 EeV. 

In Ref. [24] , we discussed the possible different cases of UHE proton images in time, angle and energy, and 
how, in each case, qualitative information could be gained on the magnetic field and the origin of UHE 
CRs. Here, we will discuss how each physical parameter can be reconstructed, and in which case. The 
physical parameters that govern the UHE CRs images are: the time delay Tioo, normalized at lOOEeV, 
the coherence length l c , the power spectrum index n B , the distance D, the emission timescale T s , the 
differential injection index 7, and the fluence N . The time delay is given by [23]: 

TE ~ L4 + (sOMpc) (lOO Eev) (l0- n G) (lMp^) yI ' ^ 

Hence, information on -B rms is contained in t wo . Both the coherence length and the distance play a double 
role. The coherence length not only contributes to the time delay, it also influences the scatter around the 
mean of the t wo — E correlation [23]: if D8 E /l c <C 1, all UHE CRs have experienced the same magnetic 
field structure during their propagation, hence the scatter is expected to be very small in the absence of 
pion production; inversely, if D6 E /l c ^> 1, the scatter is expected to be significant, At e /t e ~ 60%, even 
for negligible energy loss. The distance also enters the time delay, and it also governs the amplitude of pion 
production, hence the high energy part of the spectrum. 
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A cluster is seen on the detector as a tri-dimensional image in angle, time and energy. As the Monte-Carlo 
likelihood calculation is very time and memory intensive, we only focus on the time-energy images in the 
following. Obviously, information is also contained in the angular image itself of the cluster. For instance, 
in the limit where D9e/I c <C 1, one expects to detect a single image, albeit shifted by a sytematic offset 
9e from the true location of the source, where 9 E is tied to the time delay through te — D9%/A. Below the 
GZK cut-off, its angular size A9/9 <C 1. Note that, provided the cluster is seen at different energies, and 
9e is greater than the angular resolution, the zero-point for 9e can be reconstructed, as 9e oc E^ 1 . In the 
opposite limit, D9e/I c 1) one expects the image to be centered on the source, with an r.m.s. angular 
size 9 E - In the intermediate limit, one expects to detect several images. Moreover, if 9 E can be measured, 
it provides an estimate of the combination DB 2 m J c . 

Main features of the time-energy images of clusters of UHE protons are described in detail in Ref.[24]. 
We summarize these results briefly, as they are important to the following. If both Ts < t wo , and r wo 
is small compared to T obs , arrival time and energy are correlated according to r E oc E~ 2 ; see Fig. 1. A 
source, such that t 100 T s and t wo T obs , can be seen only in a limited range of energies, at a given 
time, as discussed in Ref.[23], as shown in Fig. 2,3. Below the GZK cut-off, the width of this stripe, in the 
time-energy plane and within the observational window of length T b s , is then tied to the ratio D9 E /l c , 
as discussed above. At the other extreme, a source emitting continuously at all energies of interest here, 
i.e. with Ts r 30 and Ts T Q b s , yields a time-energy image in which the distribution of arrival time vs. 
energy is uniform, i. e. events of any energy can be recorded at any time, as shown in Fig. 4. Finally, for a 
source, such that t 100 < T s and r 30 > T s , together with T s T obs , there exists an energy Eq, such that 
te c — Tg. In this case, protons with an energy lower than Eq are not detected, as they could not have 
reached us within T Q b s , even if they were among the first emitted. However, protons with an energy higher 
than Eq are detected as for a continuously emitting source, i. e. with a uniform distribution of arrival times 
vs. energy, see Fig. 5. 

Typical simulated clusters corresponding to these five main situations are shown in Figs. 1-5. The fluence 
N was normalized in each case so that ~ 40 events are expected within 5 years. 



3 Maximum likelihood reconstruction 

In this section, we discuss, in turn, how each parameter can be obtained from a likelihood study of UHE 
CR clusters. Certain marginalizations of Eq. (1) are used whenever the focus is only on one or a part of the 
parameters. The other parameters are then averaged or integrated over, applying weight functions (i.e., 
Bayesian priors) that represent the prior knowledge on their values. As we have currently no information 
on the fluence, the emission timescale Ts and the time delay r E , the prior chosen would be uniform in the 
logarithm of these parameters. However, we note that the time delay te is bounded from above by the 
Faraday rotation data bound on B rras l l J 2 as combined with Eq. (2). Moreover, information contained in 
the angular image should also be included in the prior on t e , as r E ^ D9 2 E /A. The marginalization over 
the injection spectral index 7 is achieved through averaging with equal weights. 

Although we focus on only one source, future large-scale experiments are expected to detect a large number 
of individual sources. Obviously, this would considerably increase the sensitivity to the physical parameters. 
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3. 1 Time delay te 

Here we assume that the source is a burst, i.e. T s <C 1 yr; we will discuss the case where T s 3> 1 yr in the 
section concerning T s . 

If the time delay is small compared to the length of the observational window, the time-energy correlation 
is scanned through, and, as Fig. 1 reveals, as simple fit of r E oc E~ 2 would allow to determine the zero- 
point of emission, hence the time delay. This constitutes a measurement of the combination DB rms l^ 2 . 
Our likelihood simulations confirm that, for the cluster shown in Fig. 1 for instance, Te is obtained within 
a factor 2. The source is found to be a burst with a high level of confidence. 

When the time delay gets significantly larger than T obs , its actual value cannot be reconstructed from 
the maximum of the likelihood. This case corresponds to the clusters shown in Figs. 2 and 3. Indeed, 
the likelihood is degenerate in the parameters N and r 10 o, as it depends mainly on the rate of detection 
Nq/Atiqq, where Arioo is the scatter in time around the mean of the te — E correlation. As long as N 
is unknown, only a lower limit to Tioo, typically r wo > T b s , can be placed. The likelihood, as calculated 
for the cluster shown in Fig. 2, and marginalized over all parameters except r wo and 7, is shown in Fig. 6 
in order to illustrate this point. The distance and the coherence length cannot be readily obtained in this 
case, as we will discuss further below. Although only a lower limit could be placed on the time delay, we 
note that, when combined with the Faraday rotation bound, this would still allow to bracket the strength 
of the extra-galactic magnetic field, within less than a few orders of magnitude. 

At this point, the information contained in the angular image of the source becomes important. If the 
angular image is not resolved, this translates into an upper limit on t e /D, which may supersede the 
Faraday rotation bound, see Eq. (2), and Eq. (3) below. At the other extreme, for a sufficiently large time 
delay, 9 E should in principle be measurable, as 



Obviously, resolving the angular image would change the prior for D and r 100 ; it would sharpen the 
maximum likelihood reconstruction, notably with respect to the various scenarios discussed in Section 2. 
We have not included this angular effect in a systematic way; a quantitative treatment of the angular images 
will be the subject of a future study. We note that the angular resolutions of future UHE CR experiments 
are fractions of a degree, hence the information contained in the angular image becomes significant for 
Tioo k, 10 yr. 

3.2 Distance D 

As mentioned above, the distance enters the likelihood mainly through the amplitude of pion production. 
As long as the high energy tail of the spectrum, i.e. E > 50 EeV, can be observed, the distance is thus 
obtained with a reasonably good accuracy from the likelihood, as marginalized over Ts, Tioo, ^V"o, and 7. In 




(3) 



6 



particular, the likelihood is sensitive to the distance if the source has a large emission timescale, Ts ^> Tioo- 
The standard error is then roughly a factor ~ 2. For example, the cluster of Fig. 4 shows a factor ~ 5 
difference in the marginalized likelihood for 50 Mpc (the true value) and 30 Mpc; such a factor is a typical 
value. The difference between 50 Mpc and 5 Mpc is typically a factor ~ 15. 

If Ts <C Tioo ^ T'obs, and the range of energies seen by the detector is above the GZK cut-off, the distance 
can still be evaluated, albeit with a somewhat larger error. Typical differences in the likelihood between 
50 and 30 Mpc and 50 and 5 Mpc are factors ~ 2 and ~ 6, respectively. 

In the intermediate case where Tioo and Ts are comparable, so that r^ c = Ts for an Eq in the observable 
energy range, the sensitivity to D is the better the lower Eq, albeit not very strong. The difference in 
the likelihood between 50 Mpc and 30 Mpc is typically a factor 3 or less (e.g., for the cluster shown in 
Fig. 5). It quickly rises, however, to a factor ~ 20 for clusters of the order of 100 events. We note that, 
due to the comparatively limited energy range seen in this case, there is a partial degeneracy between D 
and the injection spectrum parametrized by 7. For example, the marginalized likelihood does not change 
significantly if D is decreased and 7 is increased (i.e. a softer injection spectrum is assumed) at the same 
time. 

Other cases, e.g., as shown in Figs. 2 and 3, do not allow to reconstruct D. 
3.3 Emission timescale Ts 

If the emission timescale is larger than the width of the observational window, the likelihood becomes 
degenerate in the ratio N /Ts, and only a lower limit to Ts can be obtained, typically Ts > T b s - However, 
if the time delay at some intermediate energy, between say 30EeV and 100 EeV, is sufficiently large, 
and comparable to the emission timescale, then both the time delay, and the emission timescale, can be 
measured as long as a lower energy cut-off is visible above which the emission appears continuous. This 
case corresponds to the cluster shown in Fig. 5. Indeed, if the time delay is sufficiently large, then 9e can 
be observationally measured according to Eq. (3). As discussed above, the likelihood has some sensitivity 
to the distance as long as events are observed over a reasonable range of energies. This sensitivity also 
depends strongly on the statistics of the UHE CR cluster. Since t e ^ D9 E /A, t e is obtained, and, as 
discussed in Ref. [24], the emission timescale corresponds to the time delay at the cut-off energy E c , below 
which no UHE CR are recorded within T Q ^ S , as follows from the definition of this lower cut-off energy, 
te c = T s . 

In reality, however, the image observed in such a situation will appear as a burst with a large time delay most 
of the time: for E < Eq, the image is similar to that of a burst with a large time delay, as t e > Ts ^> T Q b s , 
i.e. only a limited range in energies is detected. Because r 30 > r Ec , most sources are seen at E < E c rather 
than at E > E c , where the image is similar to that of a continuous source. Notably, the likelihood for a 
bursting source with T s ~ does not exclude the above intermediate situation, for a cut-off Eq above the 
observed stripe in the time-energy image of a bursting source. In the case of Fig. 7, corresponding to the 
cluster shown in Fig. 3, the stripe is observed between ~ 30 EeV and ~ 80 EeV, and the likelihood does 
not exclude the above intermediate case with Eq > 80 EeV. Needless to say, the best reconstruction of r 100 
and T s takes place when the source is observed above E c , see Fig. 8. 
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Finally, note that if T$ ^> T30, the continuous source is hardly mistaken for a burst with a large time delay, 
which would be the closest approximation to the time-energy image of a continuous source. This can be 
seen in Fig. 9, which represents contours of the likelihood in the r 10 o — T$ plane. If the likelihood is further 
marginalized with respect to T s or r wo (see Fig. 10), a burst with a large time delay is ruled out to about 
95% confidence level. Qualitatively speaking, the difference is that for a burst with a large time delay, the 
maximum fluence occurs at some intermediate energy, and the fluence decreases with decreasing energy 
below. For a continuous source, in contrast, the fluence increases with decreasing energy, according to the 
injection negative power law spectrum. 



3.4 Injection spectrum index 7 

The injection spectrum index 7 can be measured provided UHE CRs are recorded over a bandpass in 
energy that is sufficiently broad. More precisely, in the case of a continuous source, i.e. T s ^> r 30 , 7 can be 
measured with an absolute accuracy of ~ 0.3. This is based on the likelihood as marginalized over T s , r 100 , 
and iVo, albeit for a known distance D. For example, for a continuously emitting source at D — 50Mpc 
with 7 = 2.0 (see, e.g., Fig. 4), we obtained a difference in the likelihood for 7 = 1.5 and 2.5 of a factor of 
about 30 and 2, respectively, on average. An example for this situation is given in Fig. 10. In the case of a 
continuous source with a time delay comparable to the emission timescale, i.e. such as shown in Fig. 5, the 
respective factors are about 2 and 1, and therefore hardly significant. For a burst with a small time delay 
such as in Fig. 1 these factors are about 10 and 1. A burst with r 100 T obs in which case the signal would 
be spread over a large range in energy, is even less sensitive to 7. In general, therefore, it is comparably easy 
to rule out a hard injection spectrum if the actual 7 > 2.0, but it is much harder to distinguish between 
7 = 2.0 and 2.5. 

Our analysis of the sensitivity to 7 was restricted to a fixed distance D, mainly because of CPU time 
limitations of the present serial version of our code. We expect that in the absence of information on D, an 
additional marginalization over D would decrease the sensitivity to 7. In particular, we already mentioned 
at the end of Section 3.2, a degeneracy of the likelihood between 7 and D for the intermediate case, where 
Ts and te c are comparable. 



3. 5 Fluence N 

Because of the degeneracy of the likelihood in N /T s and/or N / At wo for large timescales, it is in general 
not possible to reconstruct N . A possible exception is the case where all the particles are detected, i.e. 
Two ^ Tjbs, and T s < T obs . 

We take advantage of this section to detail slightly the marginalization procedure over N . In most cases 
we marginalized over the fluence analytically, noting that the dependence of the likelihood on N can be 
written as 

In C = a exp(x) + bx + c , (4) 
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with x = In A^o and where a, 6, and c depend on all other parameters except Nq. This just follows from the 
fact that pj in Eq. (1) is proportional to N . By using the approximation 



\n£(x) ~ ln£ 



2 V*^ -^maxj 



(5) 



in terms of value £ max and location a; max of the N — maximized likelihood, marginalizing over A^ with a 
uniform prior for In N then amounts to computing 



3. 6 Coherence length l c 

Our simulations confirm the suggestion of Ref.[23], that the main effect of l c on the angle-time-energy image 
comes through the relative size of the scatter around the 9e — t~e — E correlations. For a time-energy image, 
if the source is continuous, with Ts ^> r 30 , the correlation between te and E is drowned in the uniform 
emission of particles within the timescale Ts, and l c plays no role. The coherence length can therefore be 
estimated only when r 100 T s , and r wo ^> T obs . In this case, the signal corresponds to that shown in 
Figs. 2 and 3, and the width of the signal is related to the value of D9e/I c - The likelihood marginalized over 
Nq and 7 for the cluster Fig. 3, assuming two different coherence lengths, l c ~ 0.25 Mpc and l c ~ 1 Mpc, is 
shown in Figs. 7 and 11. The qualitative behavior of these contour plots is the following. The actual value 
of l c used in simulating the cluster is l c ~ 0.25 Mpc, so that the likelihood shown in Fig. 7 uses the correct 
value of l c . The width of a stripe in the time-energy plance, is tied to the ratio D9e/I c , or, equivalently, to 
(Dt e ) 1/2 /Ic- The likelihood shown in Fig. 11, assuming l c ~ 1 Mpc, is thus similar to that corresponding 
to l c ~ 0.25 Mpc, although it requires a comparatively larger (i.e. ~ (1/0. 25) 2 times larger) time delay to 
reproduce the large scatter in the stripe. 

This demonstrates that for a broad observed energy dispersion, a large coherence length can be ruled out 
at least when some information on the distance D and the deflection angle 9e, and thus on te, is available. 
In contrast, ruling out a small coherence length for a small observed energy dispersion is much harder, due 
to the nature of Poisson statistics. 

As mentioned previously, provided te ^ 10 4 yr, 9e can be directly measured. Note that the upper 
limit on the magnetic field strength, obtained through Faraday rotation measurements, implies te ^ 
2 x 10 5 yr(D/10 Mpc) 2 (i?/10 EeV) -2 . Whenever 9e is measurable, the degeneracy of the likelihood, with 
respect to the width of the signal in energy, thus concerns only the ratio D/l c . As we argued above, the 
likelihood itself is not sensitive to the distance if the mean energy of the signal lies below the GZK cut-off. 
If the mean energy lies above the GZK cut-off, then not only 9 E should not be measurable, but the width 
itself of the signal now arises from two very different effects: one due to the different trajectories followed by 
different UHE CRs through the magnetic field, another due to the pion production stochastic broadening 
of the signal; these effects cannot be easily disentangled. Hence, the only way this degeneracy between D 
and l c could be broken is through the observation of an astrophysical counterpart to the source, or the 
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source host, and the direct measurement of its distance. 

Finally, the likelihood was found extremely insensitive to the index of the power spectrum of magnetic 
inhomogeneities. 



4 Conclusions 

We have presented a Monte-Carlo likelihood analysis of the potential of future large-scale UHE CR exper- 
iments to reconstruct the physical parameters of the source and of intervening magnetic fields, when the 
strength of the latter is sufficient (-B rms <^ 1CT 12 G) to affect the propagation of UHE CRs. 

We discussed five generic situations of the time-energy images of UHE CRs, which we classify according 
to the values of the time delay te induced by the magnetic field, the emission timescale of the source 
T s , as compared to the lifetime of the experiment. For each case, we simulated clusters of UHE CRs 
using the instrumental characteristics typical of future experiments such as the Telescope Array, the High 
Resolution Fly's Eye, and most notably, the Pierre Auger Project. To this end we have simulated the 
emission, the propagation in an extra-galactic magnetic field, and the detection of clusters of UHE CRs. 
We then performed a Monte-Carlo likelihood analysis on these UHE CRs, and tried to reconstruct the 
physical parameters from the maximum of the likelihood. We simulated clusters of ~ 40 events, as the 
next generation experiments are expected to detect ~ 20 — 100 events per cluster if the clustering recently 
suggested by the AGASA experiment [3] is real. 

In summary, the likelihood presents different degeneracies between different parameters, which complicates 
the analysis. As an example, the likelihood is degenerate in the ratios N /Ts, or iVo/Arioo, where N is 
the total fluence, and Arioo is the spread in arrival time: these ratios represent rates of detection. Another 
example is given by the degeneracy between the distance D and the injection energy spectrum index 7. 
Yet another is the ratio (Dte) 1 / 2 /l c , that controls the size of the scatter around the mean of the t e — E 
correlation. Therefore, in most general cases, values for the different parameters cannot be pinned down, 
and generally, only domains of validity are found. 

We find that the distance to the source is obtained from the pion production signature, above the Greisen- 
Zatsepin-Kuzmin cut-off, when the emission timescale of the source dominates over the time delay. Since 
the time delay decreases with increasing energy, we find that the lower the energy Eq, defined by te c — Ts, 
the higher the accuracy on the distance D. The error on D is, in the best case, typically a factor 2, for 
one cluster of ~ 40 events. In this case, where the emission timescale dominates over the time delay at all 
observable energies, information on the magnetic field is only contained in the angular image, which we 
did not systematically include into our likelihood analysis due to computational limits. A more detailed 
investigation of angular images will be presented separately in a forthcoming study. Qualitatively, the size 
of the angular image is proportional to B Tms (Dl c Y^ 2 , whereas the structure of the image, i.e. the number of 
separate images, is controled by the ratio D 3 ^ 2 B rras /iy 2 . Finally, the case where the time delay dominates 
over the emission timescale, with a time delay shorter than the lifetime of the experiment, also allows to 
estimate the distance with a reasonable accuracy. 
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The strength of the magnetic field can only be obtained from the time-energy image in this latter case be- 
cause the angular image will not be resolvable. When the time delay dominates over the emission timescale, 
and is, at the same time, larger than the lifetime of the experiment, only a lower limit corresponding to 
this latter timescale, can be placed on the time delay, hence on the strength of the magnetic field. When 
combined with the Faraday rotation upper limit, this would nonetheless allow to bracket the r.m.s. mag- 
netic field strength within a few orders of magnitude. Here as well, significant information is contained in 
the angular image. 

The coherence length enters the ratio (Dte) 1 ^ 2 /l c that controls the scatter around the mean of the r E — E 
correlation in the time-energy image. It can therefore be estimated from the width of this image, provided 
the emission timescale is dominated by te (otherwise the correlation would not be seen), and some prior 
information on D and te is available. 

As a concluding remark, we point out that the magnetic field, although it 'scrambles' the images of UHE 
CRs, also brings extra-information. It not only leaves a signature of its own, it may also, in the case 
where the time delay becomes comparable to the emission timescale at some intermediate energy, allow an 
evaluation of the emission timescale of the source itself. There is therefore very important information hiding 
in the angle-time-energy images of UHE CRs, which could be exploited by future large-scale experiments. 
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Fig. 1. A typical cluster of UHE CR above ~ 30EeV in the time-energy plane. This cluster was produced by a 
discrete source at D = 50Mpc with an emission time scale T§ <C lyr (i.e. a burst), an average time delay of 
^ioo = 0.3 yr, and with iVo = 40. For the extra-galactic magnetic field, a power law index ng = and a coherence 
length l c ~ IMpc was assumed. The error bars correspond to an energy resolution AE/E = 0.14. The top label 
indicates the number of events in the cluster. 
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Fig. 2. Same as for Fig. 1, except for the parameters rioo = 50 yr, and A^o = 4 x 10 3 . This serves as an example for 
a burst with a long time delay, but still in the limit of small deflection, D6e/I c 1, leading to a small detected 
energy dispersion. 
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Fig. 3. Same as Fig. 1, except for the parameters rioo = 250 yr, No = 2 x 10 4 , and l c ~ 0.25 Mpc. This serves as 
an example for a burst with a long time delay for intermediate deflection, D9e/I c ~ 1- The distinct sub-bands are 
due to multiple source images. 
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Fig. 4. Same as Fig. 1, except for the parameters tioo = 0.1 yr, Ts = 500 yr, and Nq = 4 x 10 3 . This serves as an 
example of a source that is continuously emitting at all relevant energies. 
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Fig. 5. Same as Fig. 1, except for the parameters tioo = 50 yr, T§ = 200 yr, and iVo = 6 x 10 3 . A lower cut-off in 
energy occurs at Eq ~ 50 EeV where te c — T$ . 
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T 100 ^ 

Fig. 6. The likelihood for the cluster shown in Fig. 2, marginalized over Ts and Nq, plotted versus the time delay 
Tioo, for D = 50Mpc and l c ~ IMpc (the true values). The solid line is for 7 = 1.5, the dotted for 7 = 2.0 (the 
true value), and the dashed line for 7 = 2.5. As a result, tioo > 5yr to about 90% confidence level. In this case 
there is clearly no sensitivity to 7, as expected. 
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Fig. 7. The likelihood for the cluster shown in Fig. 3, marginalized over iVo and 7, plotted in the T100 — ?s plane, 
for D = 50Mpc and l c ~ 0.25 Mpc (the true values). The maximum of the likelihood occurs for noo — a few 
hundred years, Ts ~ a few years which is a good reconstruction of the true values. The contours shown go from 
the maximum down to values of about 0.01 of the maximum in steps of a factor 10 2 = 1.58. Note that values 
in the range te = T3 with E > 80EeV and Ts > 10 yr are not significantly excluded, as expected (see text). The 
fall-off at T100 50 yr and Ts < 3yr is a numerical artifact due to the limited number of propagated particles 
(4 x 10 per parameter set) which causes too patchy histograms in arrival time. 
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Fig. 8. Same as Fig. 7, but for the cluster shown in Fig. 5, and assuming l c ~ 1 Mpc (the true value). The maximum 
of the likelihood is near rioo = 4yr, Ts = 100 yr, but any values along the ridge defined by Ts ~ T50 are roughly 
equally likely, as expected from the fact that Eq ~ 50EeV (see Fig. 5 and text). 
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Fig. 9. Same as Fig. 7, but for the cluster shown in Fig. 4, and assuming l c ~ 1 Mpc (the true value). The maximum 
of the likelihood is near noo = 0.1 yr, Ts = 10 3 yr which is a good reconstruction of the true values. Parameters in 
the range Ts < 10 yr and tioo <^ 10 yr, that would be typical for a burst with a large time delay, are excluded at 
about 95% confidence level (see text). 



21 



/ 

o y ~ 
J -68 - 



-70 h 

i i i i i i i i i i i i i i i 

0.1 1.0 10.0 100.0 1000.0 

T s [yr] 

Fig. 10. The likelihood for the cluster shown in Fig. 4, marginalized over tioo and Nq, plotted versus the emission 
timescale Ts, for D = 50Mpc and l c ~ 1 Mpc (the true values). The average over the curves shown for 7 = 1.5 
(solid line), 7 = 2.0 (dotted line; the true value), and for 7 = 2.5 (dashed line), therefore, corresponds to a 
marginalization of the likelihood shown in Fig. 9 over noo- This demonstrates even more clearly that Ts < 10 yr 
is ruled out at about 95% confidence level. 
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Fig. 11. Same as Fig. 7, but assuming a magnetic field coherence length l c ~ 1 Mpc. The maximum of the likelihood 
is near noo = 400 yr, Ts = 20 yr which again is a reasonable reconstruction of the true values. Note that if rioo were 
known to be smaller than ~ 50 yr, a coherence length as large as l c = 7 Mpc could be ruled out, but l c ~ 0.25 Mpc 
would be allowed (see Fig. 7 and discussion in the text). 
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